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Abstract 

This paper describes the Monte Carlo simulation developed specifically for the VCS 
experiments below pion threshold that have been performed at MAMI and JLab. 
This simulation generates events according to the (Bethe-Heitler + Born) cross- 
section behaviour and takes into account all relevant resolution-deteriorating effects. 
It determines the "effective" solid angle for the various experimental settings which 
are used for the precise determination of the photon electroproduction absolute 
cross section. 
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1 Introduction 



Virtual Compton Scattering (VCS) off the nucleoli N is a valuable reaction to 
study the structure of the nucleon. VCS refers to the reaction 7* + iV — > 7 + iV 7 , 
where 7* and 7 represent a virtual and a real photon, respectively. It can be 
seen as an extension of Real Compton Scattering (RCS) to photon virtuality 
Q 2 7^ 0. In this case six electromagnetic observables, called Generalized Po- 
larizabilities (GPs), enter the cross section and may be determined to gain 
valuable insight into the structure of the scatterer. In the real-photon limit, 
Q 2 = 0, two of the six independent GPs are proportional to the well-known 
polarizabilities a and (3 obtained from RCS. The concept of GPs has first been 
worked out by Arenhovel et al. [1] for nuclei and later by Guichon et al. [2] 
for the nucleon. 

VCS off the proton is studied using the p(e, e'p') r f reaction: an electron scat- 
ters off a proton and a real photon is produced. The scattered electron and 
the recoiling proton are detected in coincidence, each in a high-resolution 
magnetic spectrometer, and real-photon production events are identified by 
reconstruction of the missing mass, which is zero in this reaction. The real 
photon can be produced either by the incoming or by the outgoing electron 
(the Bethe-Heitler contribution to the reaction) or by the nucleon. The nu- 
cleon contribution contains the Born part and the non-Born part. The sum of 
the Bethe-Heitler and the Born contributions will be denoted by BH+B. The 
non-Born part contains the GPs, which are accessible through the deviation 
of the measured p(e, e'p') r y cross section from the BH+B cross section, the 
latter being perfectly calculable once the elastic form factors of the proton are 
known. 

The very first dedicated VCS experiment below pion threshold to obtain in- 
formation on the GPs took place at the Mainz Microtron MAMI (Mainz, 
Germany) at Q 2 = 0.33 (GeV/c) 2 [3]. For the kinematics of this experiment 
the contribution of the GPs to the cross section had been estimated to amount 
to 10% [3]. This means that the absolute cross section had to be measured 
very precisely. In addition, one needed very elaborated analysis methods. The 
present paper is devoted to the description of the latter, which have been 
developed further and adapted to analyse also the next VCS experiment, per- 
formed at the Thomas Jefferson National Accelerator Laboratory JLab (New- 
port News, USA) at Q 2 values of 0.9 and 1.8 (GeV/c) 2 [4]. Both experiments 
are unpolarized and they are very similar, in apparatus as well as in method. 
Most numerical examples given in this paper refer to the MAMI experiment. 

The Monte Carlo code simulates p(e, e'p')7 events comparable to those of 
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the experiments. The simulation generates realistic spectra in the physical 
variables of interest and it has been used to determine with great accuracy 
what we will call effective solid angle. This effective solid angle is defined such 
that it does not only represent the geometrical acceptance, but it also includes 
the convolution of many effects. The aim of the present paper is to explain how 
the cross-section behaviour and the various resolution-deteriorating processes 
taking place in the target and in the detection systems have been taken into 
account. In addition the calculation of the simulated luminosity is explained. 
Throughout the paper the cross section used in the simulation will be often 
called "VCS cross section" for simplicity. 

The paper is organized as follows: in section 2 the kinematics of the reaction, 
a description of the experiment and the definition of the effective solid an- 
gle are discussed. In section 3 we outline the method used to implement the 
cross-section behaviour and we define the phase space in which the events are 
generated. Section 4 is devoted to a detailed description of the implementation 
of the radiative effects. Section 5 discusses the simulation package. Section 6 
covers the determination of the simulated luminosity and section 7 the calcu- 
lation of the effective solid angle. Results are presented in section 8. Finally, 
section 9 is a brief summary of the paper. 



2 Introductory definitions 

2. 1 The kinematics of the reaction and the experimental realization 

— * 

In the process p(e, e'p')j an incoming electron with momentum k scatters off a 
proton by exchange of a virtual photon 7* with momentum q and a real photon 

— * — * 

with momentum q' is emitted. The vector k and the momentum vector of the 
outgoing electron, k', define the scattering plane. The momentum vector of the 
recoiling proton, p' and q' define the reaction plane. The vector q, which is 

— * — * 

determined as k — k', lies in both planes. The direction of the real photon in the 
CM-system of 7* and p is determined by the angle between the two photons, 
#77, cm, and the angle f between the scattering and the reaction plane as is 
shown in figure f (throughout this paper all variables in the center of mass 
have an index cm; if no index is given the variable is defined in the laboratory 
system), if is defined equal to 0° when q' lies in the scattering plane and points 

— * 

to the same side of q as k' . In the CM-system, 7 and p' move back to back. 
In the laboratory system the recoiling proton is boosted in a (narrow) cone 
around q, while the undetected 7 can have any direction. This very welcome 
feature of the VCS kinematics makes it possible to cover a large range in # 77iCm 
by detecting the proton within the moderate solid angle of a high-resolution 
spectrometer. 
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CM-frame 



LAB-frame 




Figure 1. The p(e,e'p')j reaction. On the left-hand side all variables are drawn in 
the center of mass of 7* and p. In the laboratory system the proton is boosted on 
a cone around q as shown on the right. 

In the experiment a monochromatic electron beam impinges on liquid hy- 
drogen, contained in a metal can of known geometry. Its temperature and 
pressure are constantly monitored. To prevent local overheating of the liquid 
(which would cause density fluctuations and as such luminosity errors), the 
beam position on the target is continuously moving using a "raster" system. 
The scattered electron and the recoil proton are both detected in magnetic 
spectrometers, the entrance collimators denning their angular acceptances. 
The electron spectrometer defines the virtual-photon acceptance. For each 
electron-spectrometer setting, several proton-spectrometer settings are used 
to cover the interesting part of the proton cone. As both spectrometers usu- 
ally rotate in the horizontal plane, one measures essentially around (p = 0° and 
ip = 180°; only at sufficiently high momentum transfer and low real-photon 
energy the full proton cone is covered by the acceptance of the proton spec- 
trometer. 



2.2 The solid-angle definition 

For an ideal experiment (free of resolution effects) the cross section is deter- 
mined from the number of counts detected in a given phase space bin, N exp , 
and the integrated luminosity, C exp , via 



where da/ dQ is a notation for the differential cross section and dfl represents 
an infinitesimal bin in the phase space under study. It is clear that in order to 
derive precise differential cross sections from the measured data, the solid an- 
gle of the detection apparatus has to be accurately known. Using equation (1) 




(1) 
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one determines the cross section averaged over the solid angle AQi, the latter 
one being a purely geometrical quantity. When the cross section has a curva- 
ture, ascribing the average cross section to the mean kinematics results in a 
bias. One can solve this bias by ascribing the measurement to an appropriate 
different kinematics (c.f. [5]). This is, however, unpractical in our case because 
the cross section depends on five kinematical variables (see section 3). In this 
case one can stick to the central kinematics (or choose any other kinematics 
in the bin) and apply an appropriate correction to the average cross section 
in order to get an unbiased result. We choose to include this correction factor 
in the solid angle by defining another solid angle Af2 2 : 




where (da/ df2)o is the cross section at the chosen point. The solid angle 
Af2 2 deviates from AQi by the amount uj, which depends on the curvature 
of the cross section over the bin and the chosen point in the bin. To obtain 
Af2 2 one must know with sufficient accuracy the cross section behaviour of 
the process under study in the phase space region under consideration. In 
principle, this must be the cross section which one is going to measure and 
which is therefore unknown at the moment of the simulation. A sufficiently 
good approximation, however, is the BH+B cross section, since it is expected 
to deviate by less than 10% from the complete p(e, e'p')7 cross section; in 
particular its curvature, which is the decisive feature in this context, should 
be a very good approximation to the real one. 

The solid angles Af2i or Af2 2 must incorporate not only the actual detection 
geometry but also the various resolution effects. This is why these solid angles 
are called "effective" and why they can only be calculated by a Monte Carlo 
simulation. 

The present simulation is used to calculate Af2 2 of the experimental setups 
used in VCS experiments and, at the same time, to compare experimental 
and simulated data on an absolute scale. To this end, one introduces a sim- 
ulated luminosity, C sim , equivalent to the experimental one. This simulated 
luminosity is defined by 



N'- 

p sim /q\ 

where (da/ dO) stands for the differential cross section in the simulation ave- 
raged over a well-known solid angle, AQ', and N' sim is the number of events 
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generated in AQ'. Once the quantity C S i m is known, one calculates the ef- 
fective solid angle AQ (which can be AQi or AQ 2 ), i n full parallellism with 
equation (1) and (2) using 



An = ^%, (4) 



sim dn 



where da/ dQ is the cross section used in the simulation and N sim the number 
of events in AQ. Sections 3 to 6 describe how the various terms of this equation 
are obtained. 



3 Cross section behaviour and phase space definition 



3. 1 The implementation of the cross section behaviour 



As mentioned above, the calculation of AQ 2 needs as input the cross sec- 
tion behaviour. The BH+B cross section, d 5 a/ dk' dQ e r dfi 77)Cm , depends on 
the variables (k, k', 9 e r, 6> 77jCm , ip), where k, k', ... are the moduli of the corre- 
sponding three-vectors. Instead of k, k! and 9 e i one can also use q cm , q' cm and 
the photon polarisation, e, which ensures that, in the cross-section grid used 
by the simulation, only the real physical space is covered. An example of how 
the cross section behaves as a function of # 77iCm and tp for fixed q cm , q' cm and e 
for the MAMI kinematics is shown in figure 2. It is symmetric with respect to 
the scattering plane, and therefore only a "half-sphere" is shown. One clearly 
observes the two peaks corresponding to real photon emission around the in- 
coming and outgoing electron directions. Over the complete angular range the 
cross section varies by orders of magnitude, but in the phase space of interest 
(which is away from the peak region), the cross section flattens substantially. 
This allows one to choose a reasonable upper limit, or envelope value, for the 
cross section sampling, which cuts through the peaks. 

The number of events in an infinitesimal phase space bin is given by 

d^a 

dN bin = C - d/c'dcos(6» e /)dv? e /dcos(6» 77iCm )dv?, (5) 

where <p e i is the angle between the scattering plane and the horizontal plane 
containing the axis of the spectrometers. To generate counts in the phase space 
according to equation (5) one uses the acceptance-rejection method [6] in five 
dimensions with a constant as envelope for the cross section. However, the 
theoretical code [7] used to calculate the BH+B cross section is too slow to 
be used on an event-by-event basis in the simulation. To solve this problem, 
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Figure 2. The five-fold differential cross section for the p(e, e'p')j reaction as a 
function of # 77 , C ra and ip (q cm = 600 MeV/c, q' cm = 45 MeV/c and e = 0.62). 

the theoretical code has been used to calculate the BH+B cross section at the 
nodes of a five-dimensional grid in the variables (q cm ,q' crn ,e,9 yytCm ,ip). Then, 
in the simulation, the cross section value is obtained by interpolating in this 
grid, which makes the calculation faster by a factor of about 1000. In practice, 
a logarithmic interpolation is performed, reaching an accuracy of better than 
1%. 



3.2 The phase space definition 

The events have to be generated according to the five-fold differential BH+B 
cross section in a phase-space volume Ak' ■ AQ e > ■ Af2 77jCm . For an efficient 
simulation, one wants to optimize this phase space. While being not too large, 
it must cover the full acceptance of the apparatus, taking into account all 
resolution effects. The following ranges in the above mentioned variables are 
used: 

• AQ e r = Acos(8 e >) ■ A(p e r. the maximum and minimum values of 9 e > and yv 
are determined taking into account the shape of the extended target, the 
position of the spectrometer and the shape of its entrance collimator and 



7 



multiple scattering effects. 

• Ak': the lower bound is given by the lower limit of the momentum ac- 
ceptance of the electron spectrometer. The upper bound is given by the 
maximum momentum of elastically scattered electrons in the AQ e / bin de- 
fined above. This upper bound is fixed independently of the position of the 
elastic line relative to the electron spectrometer's momentum acceptance, 
since an electron, scattered with a momentum larger than the maximum 
accepted momentum, can still be detected due to energy losses before the 
spectrometer's entrance. 

• Af2 77iCm : the outgoing photon can go in any direction in the CM-system, 
therefore events are generated in the full solid angle 4n. As a result, the 
outgoing proton is also sampled in its full angular phase space, i.e. the 
full proton cone in the laboratory. This ensures that all detectable events 
are indeed taken into account, even with resolution-smearing at the target. 
Another advantage is that the simulation can be run for several proton- 
spectrometer settings all at once. For each generated proton, the simulation 
performs a loop over the various proton-spectrometer settings and tests if 
the particle is accepted or not. 

The five-fold differential cross section depends on the incoming electron mo- 
mentum at the interaction point, k, but it is not differential in this variable. 
However, although the beam is monochromatic (at the 10~ 4 level), k is not 
a constant. Each incoming electron loses energy in the target by collisions 
and by external bremsstrahlung in the material before the vertex point and 
by internal bremsstrahlung at the vertex point itself (see section 4.1). The 
resulting distribution of k at the interaction point is depicted on figure 3 for 
an incoming electron momentum of 766.4 MeV/c and a hydrogen target of 
340 mg/cm 2 . 



4 The radiation tail 

4-1 The necessity to simulate a radiative tail 

The radiative tail is a well-known feature of electron scattering experiments: 
after correction for the energy losses by collisions, the energy spectrum of the 
scattered electron shows a peak at the kinematically expected value, but this 
peak is accompanied by a radiation tail to lower energies [8]. This tail is due 
to energy loss of the incoming and outgoing electron via ionisation, external 
bremsstrahlung in the materials of the target and up to the spectrometer's 
entrances and via internal real radiation in the scattering process itself. These 
effects are of course also present in VCS experiments and give rise to the radia- 
tive tail observed in the spectrum of the missing mass squared M|, defined as 



8 



10 5 
10 4 

o 

O 1Q 3 



10 2 



100 200 300 400 500 600 700 800 

k at interaction point (MeV/c) 

Figure 3. The momentum distribution of the VCS inducing electrons as ob- 
tained by the simulation for a MAMI beam momentum of 766.4 MeV/c and 
Q 2 = 0.33 GeV 2 /c 2 . The energy losses by collision and by external and internal 
bremsstrahlung (before scattering) are taken into account. The insert shows a re- 
binned zoom on the peak region in linear scale. 

(k + p — k' — p') 2 (bold characters represent the four- vector of the particles). 
The resulting tail is shown in figure 4. 

For the calculation of the effective solid angles, one needs a recipe to generate 
in the Monte Carlo simulation the radiation tail as observed. Indeed, experi- 
mentally one applies a cut in the M\ spectrum around to select real-photon 
production events, and the same cut must be applied to the simulated events. 
The simulation reproduces the radiative tail well, which is very important 
because one wants the final cross-section result to be independent of the cut 
in M\. In fact, the influence of the position of the cut in the missing mass 
squared on the resulting cross section was lower than 1 % in the MAMI case. 
By reproducing the radiation tail in the simulation, the part of the radiative 
corrections which changes the kinematics of the reaction is taken into account, 
and the simulated radiative tail is properly convoluted with the detector ac- 
ceptance (these points will be discussed below). The other part of the radiative 
corrections is applied as a constant factor to the calculated cross section. 

Internal and external real radiation are incorporated in the simulation. These 
processes are simulated by sampling in an energy-loss distribution for the in- 
coming and outgoing electron. In the simulation only the electron's energy 
is changed, while its direction is assumed to be unaffected by the radiation 
effects (angular peaking approximation). For ionisation not only energy losses 
are taken into account for the electron and proton, but also multiple scatter- 
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Figure 4. The experimental (solid) and simulated (dashed) distributions of the 
missing mass squared, M|, for one of the MAMI kinematics (q' cm = 90 MeV/c, 
£ = 0.645, q cm = 600 MeV/c). For the simulation the BH+B cross section was used 
and the simulated distribution is normalised using the factor C exp / ' C s im- 

ing is incorporated in the simulation. The used probability distributions are 
discussed in the following subsections. 



4-2 Ionisation and multiple scattering 



Collisions of the particles in the materials of the target are simulated by apply- 
ing an energy loss and a scattering angle. The program glando of the CERN- 
libraries [9] is used to generate a realistic energy-loss distribution based on 
the mean value of the energy loss, which is calculated using the Bethe-Bloch 
equation. The deflection caused by multiple scattering is treated as explained 
in [10]. 



4-3 External radiative effects 



An electron passing through a slice of material of thickness t (in units of 
radiation length) emits photons due to bremsstrahlung. The energy loss of the 
electron, AE, is equal to the sum of the energies of all produced photons. The 
distribution of AE is given in very good approximation by [11] (t < 0.05) 



Iext(E , AE, t) 



bt 



0.57726t 



/ ae \ bt r i 



AE 



AE 3,AE., 



E 4 E 



(6) 



Eq is the kinetic energy of the electron before bremsstrahlung and b — \. 
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To generate the energy loss of the electron, one samples an energy loss ac- 
cording to the distribution (6) using the acceptance-rejection method, using 
as an envelope. To avoid variable overflows in the code for very small AE, 
the introduction of a lower limit, AEu, is necessary. In the present simulation 
AEu = 1 keV (which is well below the resolution of the experiment). Finally 
the electron energy is decreased by the obtained value for AE. 

In the simulation of the bremsstrahlung, only the energy of the electron 
is changed, which is equivalent with photon emission along the electron- 
momentum direction. This is a good approximation, since bremsstrahlung is 
very forwardly peaked. The smaller AE, the better this approximation. More- 
over, the scattering angle due to bremsstrahlung is small compared to that 
from multiple scattering. 



4-4 Internal radiative effects 

4-4-1 Virtual and Real Internal corrections 

The cross section for the p(e, e'p')^ reaction, ath, i.e. a process involving only 
one virtual photon and one real photon, cannot be measured directly, since 
in reality the pure p(e, e'p')7 process is always accompanied by additional 
photons, either real or virtual. These internal radiative effects give rise to a 
measured cross section, a exp , which deviates from a t h'- 



The correction term 5 tot is negative and depends on the cut in the radiative 
tail accompanying the scattering process. The internal radiative corrections 
to VCS are discussed in great detail in [12]. Written in first order, one gets 



5 vac accounts for vacuum polarisation diagrams, 5 ver is the vertex correction 
and S ra d is the correction for radiation in the one additional photon approxi- 
mation. One can approximately take into account higher order radiative cor- 
rections by writing [12]: 




(7) 




vac 



~\~ Over r ad i 



(8) 



0", 




(9) 



For Q 2 » m 2 , one can write: 
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a(, f(AE 



'rad ' 



7T 



In 



PC ^ 

cm) 



E E 

^cm-L^cm 



In 



\m? 



1, »/E, 



--In 



2 / J - J cm 



2 



*--£H+KS)}- 

where _E cm (^m) i s ^ ne incoming (outgoing) electron (kinetic) energy at the 
reaction vertex, a is the fine-structure constant and m is the electron mass. 
Sp is the Spence function, eg. [12]. The virtual correction terms 5 ver and 5 vac 
are independent of the cut in the radiative tail, AE c cm , and nearly constant 
for the phase space of interest. The correction for these effects will be applied 
by a constant correction factor to the measured cross section. Since only the 
first term of 5 ra d is dependent on AE^ m , this term is related to the radiative 
tail. The other terms of 5 ra d are independent of the cut position and they can 
be considered to be constant over the phase space of interest. Therefore they 
will be treated in the same way as 5 ver and 5 



' vac. ■ 



The radiative tail appears in the spectrum of the missing mass squared M\. 
The cut position should be expressed in terms of M\ since in the experiment 
one cuts in M\ to identify photon-production events. The relation between 
AE c cm and M\ is given by [12] 



J M x 

A£L = (13) 

Given the relationship (13) one could apply the correction (9) to obtain a t h, 
without including the internal radiative effects in the simulation. This proce- 
dure would only be valid if the acceptance of the detectors would not cut in 
some parts of the phase space more severely in M\ than the cut on the miss- 
ing mass itself. This, however, is not the case in the experiments. Therefore 
the simulation must generate the full radiative tail by implementing electron 
energy losses by radiation, reproducing in this way realistic spectra. 



4-4-2 Generating a radiative tail due to internal real radiation 

The first factor of the correction factor e 6rad is the product of a number of 
factors, of which the first one can be written as 

f (AgyV = /AE^y/AE^y 
V E cm E' cm J V E cm J V E' cm J 
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where a — - 

TV 



ln^^) — 1 • Assuming angular peaking, we can write [12] 
AE c cm y ( AE c cm \« /ABAVAF; 



E \ E' J \ E \ E' 



(15) 



Following [12] we interpret the factors (AE e /E e ) a and (AE' e /E' e ) a as the frac- 
tion of incoming and outgoing electrons respectively, which have lost less than 
AE e or AE' e due to internal real radiation. To sample each of these energy 
losses AE one uses the distribution, I int (E, AE, a), such that: 



I 



Integration yields 



AE / AE\ a 

I mt (E, AE, a) d(AE) = ^— J . (16) 



/ iBt (£,A£ ) a) = ^(y] , (17) 



which is normalised to 1. Remark the similarity between Ii nt (E, AE, a) and the 
leading term of I ext (E, AE, t) (eq. (6)). bt has been replaced by the quantity a, 
which is well known in literature as equivalent radiator [13], i.e. an imaginary 
radiator placed before and after the scattering center to generate internal real 
radiation. 

The recipe used to introduce the radiation tail due to internal radiation in the 
Monte Carlo simulation is: 

(1) Sample an energy loss, AE e , according to the distribution (17) with E = 
incoming electron energy E e . 

(2) Generate the kinematics of a p(e, e'p')7 event at the vertex (see figure 1) 
for the reduced energy E e — AE e of the incoming electron. The events are 
sampled according to the cross section at this reduced energy. After the 
scattering process the outgoing electron has an energy E' e at the vertex. 

(3) Sample an energy loss, AE' e , according to the distribution (17) with E = 
E' e . The outgoing electron energy is now E' e — AE' e . 

Remark that the above procedure implies electron-energy losses both at the 
incoming and the outgoing electron sides, which is fully consistent with the 
exponentiation idea. Again, for numerical reasons, one has to introduce a 
Ai?; r value to sample in the I int distribution. In practice the sampling is done 
uniformly in the integrated distribution of I int , then solving analytically for 
AE. To calculate the equivalent-radiator thickness a, one needs the value of 
Q 2 for the event, which one can only calculate after the complete process has 
taken place. However, due to the slow variation of ln(^), one gets a very 
good approximation by using the value of Q 2 given by elastic electron-proton 
scattering at the nominal beam momentum ki and scattering angle 9 e . 
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5 The simulation package 



The simulation consists of three separate programs: vcssim, resolution and 
analysis. The first one, vcssim, generates p(e, e'p')^ events in the target, 
applying ionisation energy losses and multiple scattering to all charged parti- 
cles and radiative effects to the electrons. The outgoing electron and proton 
are tracked up to the entrance of the spectrometers, where the collimator- 
acceptance cut is applied. This program produces two output files: one con- 
tains the generated events and the other one contains statistical information. 
The second program, resolution, applies the resolution effects of the spec- 
trometers on the events generated by vcssim, producing a datafile with the 
events affected by the spectrometer resolution. The third program, analysis, 
analyses the output datafile from resolution in the same way experimental 
data are analysed and produces a third datafile. The latter contains a set of 
reconstructed variables to be compared to the experimental ones. The mod- 
ular structure of the package has the advantage that one can change e.g. the 
spectrometer-resolution effects or the analysis, without having to redo the first 
step, which is the most time-consuming one. The three programs are described 
in more detail below. 



5.1 Vcssim 

Using all necessary input parameters, the program first defines the phase space 
in which it is going to sample (see section 3.2). In order to obtain an event the 
following steps are taken: first the transverse beam position on the target is 
generated in a horizontal and vertical distribution similar to the experimental 
one. An interaction point along the beamline is chosen uniformly inside the 
target length. The incoming electron is subject to multiple scattering, energy 
loss by collision and external bremsstrahlung in the target wall and the liquid 
hydrogen till the reaction vertex. Then the real internal radiation at the VCS 
vertex is simulated by an additional energy loss of the incoming electron using 
the equivalent-radiator approach discussed in section 4.4.2. Then the four- 
vector k of the electron inducing the actual VCS process is obtained. The 
energy loss through radiation can be so large, that k can already be too small 
to give any detectable electron in the final state. At this fixed value of k, the 
highest value of k! is given by the kinematics of the elastic process ep — > e'p'. 
So at this point a test is made if the momentum of the elastically scattered 
electron is high enough to be accepted in the electron spectrometer. If the test 
is negative, the event is terminated, and a new event is generated starting all 
over again. 

If the test is positive, one generates a scattered electron in the labframe and 
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an outgoing real photon direction in the CM frame. The variables cos# e /, </v, 
cos0 77;Cm , tp and k! are all sampled uniformly in their phase space. Remark 
that the outgoing real photon energy is already determined by the electron 
kinematics. If the generated kinematics is physically possible, the cross section 
is calculated for this event by interpolation in the BH+B grid of section 3.1. 
With this value for the cross section one samples a random number between 
zero and the envelope value. If the value is higher than the calculated cross 
section, the event did not pass the acceptance-rejection test and the event is 
terminated. 



As a next step, one has to determine whether the scattered electron and out- 
going proton enter the acceptances of the spectrometers. To this end, the 
momenta and directions of the electron and proton have to be calculated. 
Based on the variables 6 e >, (p e > and k', one can immediately calculate the four- 
vectors of the scattered electron and the virtual photon. Then the momentum 
four-vector in the center of mass for the outgoing real photon can be calcu- 
lated using # 77)Cm and ip. The real photon is transformed to the lab to obtain 
the four-vector q'. The four-vector of the outgoing proton can now be calcu- 
lated as p' = p + k — k' — q'. The scattered electron is first subject to real 
internal radiation energy loss. Then it loses energy by collision and by exter- 
nal radiation and undergoes multiple scattering in the various materials from 
the vertex point to the entrance collimator of the spectrometer. Similarly, the 
outgoing proton will undergo energy loss by collision and multiple scattering 
on its way to the collimator (bremsstrahlung is negligible for such a heavy 
particle). Several options are proposed to calculate the collisional energy loss 
of particles: the mean energy loss, the most probable energy loss, or a realis- 
tic energy-loss distribution (Landau distribution). For the calculation of the 
effective solid angle this last option was chosen. The spectrometer acceptance 
is defined in different ways depending on the experiment. In the case of the 
MAMI experiment the angular acceptance is defined by the collimators at the 
entrance of the spectrometers, in the case of the JLab experiment it is defined 
by cuts in a five-dimensional phase space. 



The output of the vcssim program is twofold: first a file is produced containing 
the events accepted by both spectrometers. One stores the kinematics at the 
vertex, the coordinates of the interaction point, and the momenta and angles 
of the particles at the spectrometer entrances. Also a proton-spectrometer 
index is stored, since several proton-spectrometer settings can be defined and 
filled simultaneously in one simulation run. The second output file contains 
the simulated luminosity C S i m (see section 6) and some statistical information 
regarding the simulation run. 
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5.2 Resolution 



The second program, resolution, introduces the resolution effects of the spec- 
trometers. In the experiment, for each particle seen in the set of two double 
vertical drift chambers (VDCs) the trajectory, measured in the focal plane, is 
traced back to the target using the spectrometer optics. This yields four inde- 
pendent variables at the target (the momentum modulus, two projected angles 
and one position coordinate). The accuracy obtained on these target variables 
reflects the resolution of the apparatus. The program resolution starts from 
the initial target variables (delivered by vcssim) and modifies them by adding 
the errors. Three options to realize this are discussed below. 

As a first option one can simply make use of Gaussian-distributed resolution 
effects on each target variable independently, ignoring error correlations. In 
this scheme, the difference between the initial target variable and the modified 
one is sampled in a Gaussian distribution of fixed width. 

In the experiment, the resolution effects of the VDCs will cause correlations 
in the resolution effects on the reconstructed target variables due to the spec- 
trometer optics. The second option reproduces these correlations. The consis- 
tency of the drift times with a straight line is used as estimate for the error on 
a track-by-track basis. In this way, also effects from multiple scattering within 
the chamber and from the used algorithm are effectively included [14]. After 
adding quadratically the contribution of multiple scattering in the spectrome- 
ter exit window, one obtains the total error on the detector coordinates, which 
is propagated through the known spectrometer optics to yield the error on the 
target variables. From the experimental data one can fill a four-dimensional 
histogram for each spectrometer, where each dimension corresponds to the er- 
ror on a given target variable. In this way one keeps track of error correlations 
(signs excluded) between the four target variables. The binning is chosen with 
equal width on the logarithm of the errors, which describes the distribution 
very precisely around the most probable value and sufficiently precisely in the 
long tails of the distribution, extending over four orders of magnitude, rela- 
tively to the width of the central peak. For each event the simulation samples 
in the four-dimensional histogram, yielding the width of the Gaussian error 
distribution on each target variable. Then one samples for this event in the 
obtained Gaussian distributions and one gets the modified target variables. 
This method has been applied in the analysis of the MAMI VCS experiment. 

As a third option, one can implement the resolution effects in the simulation 
directly at the detector level. In this scheme, the accepted particle is trans- 
ported to the focal plane of the spectrometer, where two types of errors are 
generated: 1) multiple scattering through the various materials, 2) the global 
resolution of the drift chambers (as deduced from experimental studies). For 
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each particle, two tracks are considered: one with and one without these focal 
plane resolution effects. As in the second option, one uses the full spectrom- 
eter optics to transport the tracks back to the target. Now the quantity of 
interest is just the difference between the two tracks for the same particle. 
This difference represents the resolution effects on the target variables. Since 
one uses the difference between the two tracks, one can approximate the op- 
tical transport from target to focal plane: it does not need to be the exact 
reverse of the optical transport from focal plane to target. The method gen- 
erates error correlations at the target, signs included. Large resolution tails 
are introduced at the level of the detector coordinates, e.g. by sampling in 
the sum of two Gaussian distributions with very different widths for the drift 
chamber resolution. This method has been applied in the analysis of the JLab 
VCS experiment. 

The output is a datafile containing the same variables as the one from vcssim, 
but now they include also the spectrometers' resolution effects. This datafile 
is comparable to the experimental one. 

5.3 Analysis 

The third and final part of the simulation, the analysis program, performs 
the full event reconstruction as in the analysis of the experimental data. From 
the reconstructed target coordinates, one first calculates the vertex point and 
from this the pathlengths of the particles in target materials and the corre- 
sponding (mean collisional) energy losses. The particle momenta are corrected 
for these energy losses, yielding the four- vectors at the vertex point. Then the 
complete reaction kinematics is reconstructed, including the missing particle. 
Then one can compare e.g. the distribution in missing mass squared M\ to 
the experimental one, as shown on figure 4. 



6 The determination of the simulated luminosity C sim 

As it is clear from equation (4) one needs to know the simulated luminosity 
in order to obtain the effective solid angle. In the experiment, the luminosity, 
C exp , is obtained as the product of the number of incoming electrons and the 
target thickness and is totally independent of the reaction under study. 

The simulation uses a different approach: the luminosity in the simulation, 
Csim, is calculated from the cross-section samples of the acceptance-rejection 
method, i.e. from the reaction itself. The method is most efficient and gives 
a very accurate result, provided the procedure is established carefully. One 
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counts the number of samples, N, generated during a simulation run in the 
luminosity phase space, L.P.S., which is a sub-part of the total simulation 
phase space (see section 3.2). Simultaneously, the cross section is integrated 
over this luminosity phase space. According to equation (3) C sim is then simply 
given by 

N 

~ 1 H », HO dB ; dfc'd^da^ (18) 

J L.P.S. dk' dQ e r df2 77jCm e 77, cm 

As such, C S im is actually calculated in a reverse way, i.e. at the end of a 
simulation run, once the number of generated events is known. 

In principle one is free to define the size of the luminosity phase space. However 
one will have to choose a luminosity phase space that is smaller than the 
simulation phase space. 

The first complication is due to the method used to implement the cross- 
section behaviour. As mentioned in section 3.1, the acceptance- rejection method 
with constant envelope is used, with a rejection level of about 90%. However, 
among these rejected samples a large fraction can be kept to calculate the 
value of the cross-section integral over the luminosity phase space. One just 
has to make sure that the luminosity phase space does not overlap with the 
regions of the simulation phase space where the cross section is larger than 
the envelope, since the acceptance-rejection method does not work in these 
regions. 

The second complication is connected to the fact that the real k distribution 
has a low-momentum tail and to the fact that the cross section depends on the 
value of k. If all VCS inducing electrons would have the same momentum k, one 
could immediately apply equation (18). This case is explained in subsection 
6.1. The case of the real k distribution, for which one cannot apply equation 
(18) directly, is discussed in subsection 6.2. 



6.1 The definition of the luminosity phase space for the case of constant k 



If all interacting electrons would have the same momentum, ki (the nominal 
beam momentum), the cross-section integral of equation (18) would be given 
by 

f d 5 cr 

I = / ,„, n , n dA;'dfi e /dfi 77)Cm . (19) 

Jl.p.s. d/c'diZe' di2 77jCm 

One has to define the luminosity phase space as an integration range in 
dn 77iCm , dfl e r and dk': this is done using a 5-dimensional box in (k', 6 e r, 
fe', #77,cm) V 9 )) where the limits on each variable are independent of the other 
variables. For example, for the MAMI experiment the box has the following 
dimensions: 



18 



k' 



elas,min 



Elastic Line 




cut 



L.P.S. 



min 



Figure 5. The luminosity phase space shown schematically in the two dimensions k' 
and 8 e i for the case of constant k. The L.P.S. is the shaded rectangle, upper-bounded 
by Kut = (Kias,min + ^min)/^ as explained in the text. The range in 9 e i is defined 
by the acceptance of the spectrometer. The solid line is the (ep — > ep) elastic line 
at the incoming energy k. 

• In Af2 77jCm : # 77iCm varies from to 7r and tp varies from about 0.8 to 5.48 
radians. This region is chosen in order to stay away from the steep cross- 
section rise in the region of the BH+B peaks around (p = (see figure 2). 

• In AQ e r. one uses the complete solid angle in which the electron directions 
are sampled. 

• In Ak'\ for the minimum of k', the lower limit of the electron spectrometer 
acceptance, k' min , is used. For the maximum of the ^'-integration range one 
has to be careful not to cross the envelope value with the cross-section 
values in the regions in Af2 77cm and in AQ e ' denned above. Indeed, as k! 
increases at fixed fcj, one approaches the elastic kinematics and as such the 
cross section rises. To stay far enough away from the elastic kinematics, 
the maximum value of k! for the integration range is taken equal to k' cut = 
(Klas,mm + Knn)/ 2 - Th e quantity k' elas>min is the minimum momentum an 
elastically scattered electron can have in AQ e >, for an incoming electron 
momentum fcj. This is illustrated on figure 5. 

6.2 Taking into account the realistic distribution of the incoming momentum 



Due to energy losses the electron momentum at the vertex point becomes 
distributed as shown in figure 3. In this realistic case one could divide this 
distribution in small bins in k and apply equation (18) to calculate the partial 
luminosity for each bin j 



k 




(20) 
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where I a j is given by expression (19) and Nj is the number of accepted sam- 
ples, both evaluated in a luminosity phase space similar to the one of section 
6.1, where the range in k is limited to the bin j. Then C sim would be equal to 
There are two limitations to this procedure: 

(1) It is not possible to use formula (20) on the whole incoming electron 
momentum distribution, because the cross-section grid does not contain 
cross-section samples for the entire — > ki incoming momentum range. 
One has to cut somewhere in the /c-range, hereby defining a cut value 

k C ut- 

(2) For small values of k, the elastic line drawn on figure 5 lies totally below 
the lower bound of the electron- momentum acceptance k' min , therefore 
k' cut lies below k' min and the luminosity phase space can not be defined as 
in section 6.1 for these electrons. 

To solve these problems the cross-section integration will be performed in a 
limited range of incoming electron momentum, i.e. a bin in k for which the 
L.P.S. can be defined as in section 6.1. This will yield a partial luminosity. The 
total luminosity will then be obtained by a simple renormalization procedure. 



6.2.1 The k-range for the cross section integration 

Quite obviously, the cross section integration of equation (19) should be per- 
formed for the incoming electron momenta that are closest to the beam mo- 
mentum ki. Therefore one defines a range of the type [k cut , ki]. The luminosity 
phase space is then defined as in section 6.1. The value of k' cut is calculated 
using the elastic line at the lowest incoming momentum of the bin, i.e. at 
k = k cu t- 

Of course, when one lowers the value of k cut , one reduces the size of the lumi- 
nosity phase space (due to the choice of k' cut ), and the statistical error on the 
luminosity increases. So one should keep k cu t close enough to ki. For example 
for the MAMI experiment the range [ki — 3 MeV/c, ki] was chosen. It contains 
about 80% of all incoming electrons, and yields a statistical error on £ sim well 
below 1 %. 

During execution the number of samples in the luminosity phase space, N L P _s., 
is counted and the integral over the cross section in the L.P.S. , I a , is calculated. 
Nl.p.s. is the number of samples accepted by the acceptance-rejection method 
of section 3.1. At the end of execution, the partial luminosity C S im,zMeV/c is 
given by N L , RS ./I a - 
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6.2.2 The renormalization factor 



By the method described above we know the luminosity C S im,3MeV/c corre- 
sponding to a fraction, /, of all incoming electrons, which have a momentum 
higher than ki — 3 MeV/c. This fraction / is easily calculated in the simulation: 
one counts the total number of A;- values that have been generated, N tot , and 
the number of values that have been generated above the threshold of ki — 3 
MeV/c, N 3 MeV/c- At the end one has / = N 3 MeV/c/N tot . However, one needs to 
know the total luminosity £ S i m according to all incoming electrons. One can 
obtain the right value by correcting for the electrons one did not count in the 
calculation of C S im,3MeV/c- Since the luminosity is independent of the reaction 
under study, the total luminosity C sim is obtained by dividing C sim 3MeV/c by 
/■ 



7 Calculation of the effective solid angle 



The data in the output file from the analysis program, in combination with 
simulated luminosity C S i m from the vcssim program are used to calculate the 
effective solid angle for any given bin in the phase space, applying equation (4): 



Afi = ^ , (21) 



-sim- dfc'dO 



where d b a sim / dk' dQ e r dfi 77iCm is now the differential cross section for the 
p(e, e'p')7 reaction, used in the simulation and N sim the number of counts 
in the bin. AQ is similar as in equation (1) or (2), with now the specific 
dimension of (sr 2 - MeV/c), as can be deduced from equation 21. By applying 
energy losses for radiative effects in the simulation, a part of the radiative 
correction is automatically taken into account in AQ. If the cross section is 
taken to be a constant value over the complete phase space, equation (21) will 
yield AQi. Calculating the cross section from the data using this AQi yields 
the experimental cross section averaged over the bin. On the other hand, if 
the simulation has been performed using the BH+B cross section, the quantity 
d 5 (J S j m / dk' dQ e ' df2 77)Cm in equation (21) equals the BH+B cross-section value 
at a given phase-space point, which can be chosen anywhere, preferentially 
in the bin. This procedure will give rise to AQ 2 , comparable to equation (2). 
Applying this AQ 2 to the measured data will result in the actual cross-section 
value at a given phase-space point. In order to get a precise result, the cross 
section in the simulation must have a behaviour very close to the true cross- 
section shape. 

If necessary, one can even use an iteration procedure to improve the value of 
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Figure 6. Effective solid angles for one of the settings of the MAMI experiment as 
a function of 7T , C m {f = n, Qcm = 600 MeV/c, q' cm = 45 MeV/c and e = 0.62). 
AOi is obtained by running the simulation with a flat cross section, for Af2 2 the 
BH+B cross section is used. A$73 is the same quantity as AS1 2 , but for a simulation 
without radiative effects. The purely statistical errors are smaller than the size of 
the symbols. 

Af2 2 by implementing in the simulation a cross section of the type BH+B plus 
a polarizability effect. This procedure was tested for the two experiments. In 
the case of MAMI, the relative change of Af2 2 was smaller than 1 %, hence 
the iterations had a negligible effect on the physics observables (the GPs). 
In the case of JLab, the relative change of Af2 2 was larger, typically a few 
percent, translating into significant changes of the physics observables: after 
the first (resp. second) iteration, the GPs reached ~ 70 % (resp. 90 %) of their 
convergence value. The full convergence was obtained after the third iteration. 



8 Results of the effective solid angle calculation 



As an example, figure 6 shows the obtained effective solid angles at q' cm = 
45 MeV/c (e = 0.62 and q cm = 600 MeV/c). The phase space is denned by 
40 MeV/c < q' cm < 50 MeV/c, 158° < ip < 202°. The statistical error on the 
effective solid angle in the plateau region is about 1%. Af^ is the solid angle 
obtained by generating the events according to the BH+B cross section. A 
simulation with a constant cross section gives Af^. It turns out that for this 
setting the difference between Af2 x and Af2 2 is up to the order of 10%. It is 
also interesting to run the simulation without radiative effects, which results 
in AQ 3 . The right panel of figure 6 shows clearly that radiative effects have to 
be included in the simulation, since there is no common scaling factor between 
Af2 2 and Af2 3 for the complete phase space. 
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9 Summary 



The Monte Carlo simulation described in this paper has been developed for the 
analysis of the VCS experiment at MAMI and has been adapted afterwards 
for the analysis of the VCS experiment at JLab. It has been used to gener- 
ate realistic observable spectra, which can be compared with the measured 
ones, and to determine accurately effective solid angles which also account 
for the radiative processes accompanying the VCS reaction. The use of a five- 
dimensional cross-section grid covering the complete simulation phase space 
allows to generate events according to the Bethe-Heitler+Born cross section at 
a very acceptable rate, using the acceptance-rejection method with a constant 
envelope. External and internal radiation of real photons are implemented in a 
well-founded way by generating realistic radiative tails and convoluting these 
effects with the acceptance of the detection system. 

The simulation described above is flexible. All resolution deteriorating effects 
can independently be switched on or off and it is possible to use a constant 
cross section or the BH+B cross section to generate events. Due to the multi- 
ple proton-spectrometer option the yield in several proton-arm settings for one 
electron-spectrometer setting can be simulated in one run, while the modular- 
ity of the code gives the possibility to study spectrometer-resolution effects in 
an efficient way. Finally, the program is general enough to allow adaptation 
to many other processes, including e.g. elastic scattering and pion electropro- 
duction. 
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